using System;
using System.Reflection;
using NBody;
using NBody.InitialConditions;
using NBody.Cosmology;
using OptionParser;

[assembly: AssemblyTitle ("NBodyCube")]
[assembly: AssemblyVersion ("1.0.*")]
[assembly: AssemblyDescription ("Generate initial perturbations for a cube.")]
[assembly: AssemblyCopyright ("2005 Joseph D. MacMillan")]

public class CubeOptions : OutputGetOptions
{
    [Option ("-n", "Number of particles per side of initial cube")]
    public int N = 32;
    
    [Option ("-l", "Length of side of box (comoving)")]
    public double L = 10000.0;
    
    [Option ("-z", "initial redshift")]
    public double Z = 50.0;
    
    [Option ("-min", "Minimum wavelength")]
    public int MinWave = 0;
    
    [Option ("-max", "Maximum wavelength")]
    public int MaxWave = -1;
    
    [Option ("-smooth", "Smoothing scale for power spectrum")]
    public double Smooth = 0.0;
    
    [Option ("-s", "Random seed")]
    public int Seed = 0;
        
    [Option ("-cosmo", "Cosmology for halo.")]
    public string CosmoString = "SCDM";
}

public class NBodyCube
{    
    public static void Main(string[] args)
    {
        CubeOptions opts = new CubeOptions();
        opts.ProcessArgs(args);
        
        Model model = Model.FromString(opts.CosmoString);
        Spectrum spec = new Spectrum(model, new TransferFunction(Spectrum.TransferBBKS));
        spec.SmoothingScale = opts.Smooth;
        spec.Normalization = Normalization.Sigma8;
                
        Cube cube = new Cube(opts.L, model.Density(0.0) * opts.L * opts.L * opts.L, Cube.Origin.Corner);
                
        Perturbations perturb = new Perturbations(spec, opts.Seed);
        perturb.MinWave = opts.MinWave;
        perturb.MaxWave = opts.MaxWave;
        Displacements[] disp = perturb.CalcGaussDispFFT(opts.N, opts.L);
        
        NBodySystem s = cube.GenerateGrid(opts.N, model.Time(opts.Z), disp);
                
        s.Write(opts.OutFile);
    }
}
